#########################################################################
# Time-stamp: <05/08/11 15:59:14 ostolop>
#########################################################################
package EP::EPC::Core::Between_Group_Analysis;

use strict;

use EP::Config;
use EP::Common;
use EP::Common::General;
use EP::Common::Error ':try';

use EP::EPC::Core;
use EP::EPC::Core::Component;

use EP::Rclient;

use File::Basename;
use File::Temp;

our @ISA = qw ( EP::EPC::Core::Component );

sub process_component {
  my $self = shift;

  my ( $output, $images, $chime ) = $self->run_bga ();

  my %process_result = ( %$output, %$images );
  $process_result{chime} = { type => "chime", value => $chime, label => "Visualization" };
  $process_result{src_dataset_id} = $self->{src_dataset}->id;
  $process_result{dst_dataset_id} = $self->{src_dataset}->id;

  $self->{CD}->{process_result} = \%process_result;
  return $self->record_analysis_history ();
}

sub nan_impute_row_avg
    {
    my ($bin_file) = @_;
    my $nanimp_bin_file = $bin_file . ".nanimp";

    return $nanimp_bin_file if -s $nanimp_bin_file;

    EP::Rclient->R(<<__REVAL__);
app_root<-"$EP::Config::h->{ROOT}";
source("$EP::Config::h->{ROOT}/cgi-lib/R/ep.io.R");
source("$EP::Config::h->{ROOT}/cgi-lib/R/ep.impute_data.R");

M <- ep.readBin ( "$bin_file" );
O <- ep.nanimprowavg ( M );
ep.writeBin ( O, "$nanimp_bin_file" );
__REVAL__

    die "R did not generate imputed data!" if not -e $nanimp_bin_file;
    return $nanimp_bin_file;
    }

sub run_bga {
  my $self = shift;
  my ($src_dataset_id) = @_;

  throw EP::Common::Error ( -value => AFFY_TRANSFORM_RAW_DATA ) if $self->{src_dataset}->type eq "affy";

  #############################
  # Retrieve training data
  ############################

  my $dataset_folder_path     = $EP::Config::h->{EPNG_USRDATA_PATH} . "/" . $self->{src_dataset}->folder->location;
  my $bin_file = $dataset_folder_path . "/" . $self->{src_dataset}->filename . ".bin";
  $bin_file = nan_impute_row_avg ( $bin_file ) if $self->{query}->param ( "impute_row_avg" );


   ########################################
   # Retrieving the dataset for suppl proj

  my @pnames = $self->{query}->param;
  my @dataset_ids = map { /chk_dataset_(\d+)/ } grep { /^chk_dataset/ } @pnames;
  my @action_ids = map { /chk_action_(\d+)/ } grep { /^chk_action/ } @pnames;

  die "Can't specify more than one dataset for supplementary projection!" if scalar @dataset_ids + scalar @action_ids > 1;

  my $testdataset;
  if ( scalar @dataset_ids == 1 ) {
    $testdataset = EP::Common::DBI::Dataset->retrieve ( $dataset_ids[0] );
  } elsif ( scalar @action_ids == 1 ) {
    $testdataset = EP::Common::DBI::Action->retrieve ( $action_ids[0] )->dst_dataset;
  }
   ########################################

  my $suppl_proj = "Y" if defined $testdataset;

  my $testdataset_folder_path = $EP::Config::h->{EPNG_USRDATA_PATH} . "/" . $testdataset->folder->location if $suppl_proj;

  ####################################################
  # Class vector
  # try the long form first, and if not passed, try short form.
  # Calls sub format_classvec
   #######################################################
  my $factor = $self->{query}->param ( "factors" );
  my $matrix = $self->{src_dataset}->data->raw_data;
  my $column = $self->{src_dataset}->data->{column};
  my $nr_cols = scalar keys %{$column};

  my @cv;
  push @cv, $column->{$_}->{full_annot}->Values ( $factor ) foreach 0 .. $nr_cols - 1;
  my $classvec = join ( ",", @cv );

  die "No class vector defined. Please enter class vector in form provided." if not  $classvec;

  #######################
  # Tmpl Parameters
  #########################

  my $rand = "$$" . rand(100);

  my $tmpl_values = { html_root    => $EP::Config::h->{HTML_ROOT},
		      src          => $bin_file,
		      srcrowannot  => $dataset_folder_path . "/" . $self->{src_dataset}->filename . ".rows",
		      srccolannot  => $dataset_folder_path . "/" . $self->{src_dataset}->filename . ".columns",
		      filename     => $self->{src_dataset}->filename.".".$rand,
		      outputdir    => $dataset_folder_path,
		      classvec     => $classvec,
		      ord_type     => $self->{query}->param ( "ord_type" ),
		      outputlist   => $self->{query}->param ( "output_files" ).$self->{query}->param ( "output_graphs" ),
		      suppl_proj   => $suppl_proj . "",
		      suppl_output => $self->{query}->param ( "suppl_output" ) . "",
		      testfilename => ($suppl_proj) ? $testdataset->filename : "",
		      testsrc      => ($suppl_proj) ? $testdataset_folder_path . "/" . $testdataset->filename . ".bin" : "",
		    };

  ##################################
  # Checks before running: check that the class vector is the same 
  # length as the number of arrays (columns)
  ##################################
  # Count the number of elements in the classvector
  my $veccount;

  my @e = split /,/, $classvec;
  $veccount = scalar @e;

  # $filecount counts the number of lines in the columns labels.
  my $colcount=0;

  open(FILE, "< $tmpl_values->{srccolannot}") or die "can't open $tmpl_values->{srccolannot}: $!";
  $colcount++ while <FILE>;
  $colcount = $colcount -1;  # Exclude header line
  close FILE;

  if ($veccount != $colcount) {
    die "Please re-enter the class vector.  The length of the class vector you entered is $veccount, but there are $colcount arrays in the microarray dataset you selected. The length of the class vector should equal the number of arrays.";
  }

  EP::Rclient->R(<<__RCODE__);
app_root = "$EP::Config::h->{ROOT}";
source("$EP::Config::h->{ROOT}/cgi-lib/tmpl/R/ep.run_bga.tmpl.R");
ep.run_bga( "$tmpl_values->{src}"
          , "$tmpl_values->{testsrc}"
          , "$tmpl_values->{filename}"
          , "$tmpl_values->{testfilename}"
          , "$tmpl_values->{outputdir}"
          , "$tmpl_values->{classvec}"
          , "$tmpl_values->{ord_type}"
          , "$tmpl_values->{outputlist}"
          , "$tmpl_values->{suppl_proj}"
          , "$tmpl_values->{suppl_output}"
          );
__RCODE__

  ########################
  # Process R output 
  #########################
  my $output_www_base = "$EP::Config::h->{EPNG_USRDATA_URL}/" . $self->{src_dataset}->folder->location . "/" . $self->{src_dataset}->filename.".".$rand;
  my $output_path_base = "$EP::Config::h->{EPNG_USRDATA_PATH}/" . $self->{src_dataset}->folder->location . "/" . $self->{src_dataset}->filename.".".$rand;

  my $output = get_run_bga_output ( $tmpl_values->{outputlist}. $tmpl_values->{suppl_output}, $output_www_base, $tmpl_values->{ord_type}, $output_path_base, $classvec );
  my $images = get_run_bga_images ( $tmpl_values->{outputlist}. $tmpl_values->{suppl_output}, $output_www_base, $tmpl_values->{ord_type} );
  my $chime  = display_chime ($output_path_base, $classvec );

  return ( $output, $images, $chime );
}


sub parse_template {
  my ($tmpl_values, $src, $dst) = @_;
  local(*FILE);


  my $template = new Text::Template ( TYPE => "FILE", SOURCE => $src );
  my $text = $template->fill_in ( HASH => $tmpl_values );

  open (FILE, ">$dst");
  print FILE $text;
  close FILE;

  return 1;
}


sub format_classvec {

   # expands class vector from short form into long format
   # short form:
   # X-Y:CLASS NAME 1;Z-Y:CLASS NAME 2;Y:CLASS NAME 3

  my $self           = shift;
  my $classvec_short = shift;

  my %classes;
  my @class_defs = split /;/, $classvec_short;

  foreach (@class_defs) {

     # Check labels 1-10:A;  Split each at the :

     my ( $range, $label ) = split /:/;
     die "Couldn't format class vector. Please check the class vertor you entered. It is likely that there is a missing class label or range." if not $label;


     # Check the range 1-10 is a number

     my ($range_from, $range_to ) = split /-/, $range;
     die "First range limit ($range_from) invalid for class '$label'!" if $range_from !~ /^\d+$/;
     die "Second range limit ($range_to) invalid for class '$label'!" if defined $range_to && $range_to !~ /^\d+$/;

     push @{$classes{$label}}, [$range_from, defined $range_to ? $range_to : $range_from ];
  }

  my @classvec_array;
  foreach my $label ( keys %classes ) {
    foreach my $range_def ( @{$classes{$label}} ) {
      my ($from, $to) = ($$range_def[0] - 1, $$range_def[1] - 1);

     # Check for overlapping ranges in the class vector
     for ($from..$to) {
     if ( defined $classvec_array[$_] ) {
         die "Overlapping ranges: sample ".($_+1)." belongs to both $classvec_array[$_] and '$label'!";
     } else {
       $classvec_array[$_] = "'$label'";
     }
   }

    }
  }

 #  my $unmapped_samples = grep /^$/, @classvec_array;
 #     if ( $unmapped_samples ) {
 #     die "Incontiguous ranges defined! $unmapped_samples samples are unmapped.";
 #  }

  my $classvec_long = join (",", @classvec_array);
  print STDERR "[CLASSVEC] ", scalar @classvec_array, "$classvec_long\n" if DEBUG > 1;
  return $classvec_long;

}


sub extract_names {
  #####################################
  # Process Class vector. 
  # Not used extinct code.
  # This produces a list of all array samples names.
  # Originally intended to produce list of arraynames from which user
  # could place samples in classes using "button" list
  ####################################

  my $self = shift;
  my $user_folder = shift;

  my $datasets;
  my $count = 0;

  my $username = $self->{DDD}->get_username();
  my $EPNG_USRDATA_PATH  = $EP::Config::h->{EPDATA}."/EPCLUST_DATA/$username";

  opendir ( DATASETS, "$EPNG_USRDATA_PATH/" ) or die "Couldn't open directory $EPNG_USRDATA_PATH/ for reading";
  my (@files) = grep(!/^\./,readdir(DATASETS));
  closedir DATASETS;

  my $filename = "/data/nfs/EP/DATA/EPCLUST_DATA/guest/Train";

  open(FILE, $filename) || die "cannot open file";

  my $tag ="EXPNAME";
  my $arraynames;

  while(<FILE>) {
    if($_=~m/<$tag>(.*)?<\/$tag>.*/sig) {
      print "$1 ";
      #$arraynames = $arraynames. $1;
      $arraynames->{$count}{key} =  $1;
      $arraynames->{$count}{value} = $1;
      $arraynames->{$count}{name} = "EXPNAME";
      $count++;
    }
  }
  close(FILE);

  my $dl = &EP::EPC::Interface::DD::dynamic_list( $arraynames );
  return $dl->findnodes('/*');
}

sub extract_classes {
   #####################################
  # Process Class vector. 
  # Not used currently.
  # This produces a list of all arrays and then the user can select
  # from a list to place samples into different classes to form a classvector
  ####################################

  my $subclass = shift;
  $subclass =~ s/\s*(\-|\_|\\.\.)\s*/ .. /g;
  my @subclass = split(/;+\s*/, $subclass);
  my $outvec;

  foreach (@subclass) {
    my $i = $_;
    my @item= split(/\,\s*/, $i);
    chomp(@item);

    my $rep;

    if ($#item ==0) {
      $rep=$item[0]. " ";
    } else {
      my $class= $item[1]." ";
      my $index= $item[0];

      if ($index =~ /\.\./) {
	my @num = split(/ .. /, $index);
	my $firstnum = $num[0];
	my $endnum = $num[1];
	#print "\n start: $firstnum end: $endnum";
	$rep = $class x ($endnum-$firstnum +1);
      } else {
	$rep = $class x $index;
      }
    }
    $outvec=$outvec.$rep;
  }

  chomp($outvec);
  return $outvec;
}

sub get_run_bga_images {
 
  ##################
  # Produce output images
  #####################
  my $outlist         = shift;
  my $output_www_base = shift;
  my $ord_type        = shift;

  my %images;

 if ($outlist =~ /D/ ) {
    $images{def_graph} = { type => "image/png",
			   value => $output_www_base . "_def.png",
			   label => "Overall plot. Graph of arrays, Genes and eigenvalues of BGA" };
  }

  if ($outlist =~ /V/ ) {
    $images{eig_graph} = { type => "image/png",
			   value => $output_www_base . "_eig.png",
			   label => "Graph of eigenvalues of BGA" };
  }

  if ($outlist =~ /A/ ) {
    $images{array_graph} = { type => "image/png",
			     value => $output_www_base . "_array.png",
			     label => "Graph of arrays of BGA" };
  }

  if ($outlist =~ /G/ ) {
    $images{gene_graph} = { type => "image/png",
			    value => $output_www_base . "_genes.png",
			    label => "Graph of genes of BGA" };
  }

  if ($outlist =~ /S/ ) {
    $images{proj_graph} = { type => "image/png",
			    value => $output_www_base . "_suppl.png",
			    label => "Graph of projected supplementary co-ordinates of BGA" };
  }

  return \%images;
}

sub get_run_bga_output {

  ###########################
  #  Produce text output
  ###########################
  my $outlist         = shift;
  my $output_www_base = shift;
  my $ord_type        = shift;
  my $output_path_base = shift;
  my $classvec        = shift;

  my %c;
  $c{$_}++ foreach split /,/, $classvec;

  my $nf = -1 + scalar keys %c;

  my %output;

  $output{logfile} = { type => "text/plain",
		       value => $output_www_base . "_log.txt",
		       label => "Summary of Between Group Analysis" };

  if ($outlist =~ /E/ ) {
    $output{eigvals} = { type => "text/plain",
			 value => $output_www_base . "_eig.txt",
			 label => "Eigenvalues" };
  }

  if ( $outlist =~ /R/ ) {
    open P, $output_path_base."_ls.txt";
    open PX, ">$output_path_base"."_ls.xml" or die "CAN'T OPEN $output_path_base _ls.xml - $!";
    print PX '<?xml version="1.0"?>';
    print PX '<table>';
    print PX '<head>';
    print PX '<data>Name</data>';
    print PX "<data>Axis $_</data>" foreach 1 .. $nf;
    print PX '</head>';
    while (<P>) {
      chomp;
      my @d = split /\t/;
      print PX '<row>';
      print PX "<data>$_</data>" foreach @d;
      print PX '</row>';	
    }
    print PX '</table>';
    close P;
    close PX;
    $output{rows} = {  type  => "xml/table",
		       value => $output_www_base."_ls.xml",
		       label => "Coordinates of rows" };
  }

  if ( $outlist =~ /C/ ) {
    open P, $output_path_base."_co.txt";
    open PX, ">$output_path_base"."_co.xml";
    print PX '<?xml version="1.0"?>';
    print PX '<table>';
    print PX '<head>';
    print PX '<data>Name</data>';
    print PX "<data>Component $_</data>" foreach 1 .. $nf;
    print PX '</head>';
    while (<P>) {
      chomp;
      my @d = split /\t/;
      print PX '<row>';
      print PX "<data>$_</data>" foreach @d;
      print PX '</row>';	
    }
    print PX '</table>';
    close P;
    close PX;
    $output{cols} = { type => "xml/table",
		      value => $output_www_base."_co.xml",
		      label => "Coordinates of columns" };
  }

  if ($outlist =~ /T/ ) {
    open P, $output_path_base."_li.txt";
    open PX, ">$output_path_base"."_li.xml";
    print PX '<?xml version="1.0"?>';
    print PX '<table>';
    print PX '<head>';
    print PX '<data>Name</data>';
    print PX "<data>Axis $_</data>" foreach 1 .. $nf;
    print PX '</head>';
    while (<P>) {
      chomp;
      my @d = split /\t/;
      print PX '<row>';
      print PX "<data>$_</data>" foreach @d;
      print PX '</row>';	
    }
    print PX '</table>';
    close P;
    close PX;
    $output{cent} = { type => "xml/table",
		      value => $output_www_base."_li.xml",
		      label => "Centroids of rows" };
  }

  if ($outlist =~ /P/ ) {
    open P, $output_path_base."_suppl.txt";
    open PX, ">$output_path_base"."_res.xml";
    print PX '<?xml version="1.0"?>';
    print PX '<table>';
    print PX '<head>';
    print PX '<data>Name</data>';
    print PX "<data>Projected Axis $_</data>" foreach 1 .. $nf;
    print PX '<data>Closest Centroid</data>';
    print PX '<data>Predicted Class</data>';
    print PX '</head>';
    while (<P>) {
      chomp;
      my @d = split /\t/;
      print PX '<row>';
      print PX "<data>$_</data>" foreach @d;
      print PX '</row>';	
    }
    print PX '</table>';
    close P;
    close PX;
    $output{proj} = { type => "xml/table",
		      value => $output_www_base."_res.xml",
		      label => "Coordinates of projected supplementary data" };
  }

  return \%output;
}

sub javacols {

   ############################
   # Defines colours used by jmol. Not called.
   # ###########################

   my $numcol= shift;  #Input the number of colours required
   my @javacols;
   my @getcols;
   my @getcols1;

  ########################
  # List of colours which can be used by jmol
  #########################

   @javacols = ("aliceblue" ,"antiquewhite" ,"aqua" ,"aquamarine" ,"azure" ,"beige" ,"bisque" ,"black" ,"blanchedalmond" ,"blueviolet" ,"brown" ,"burlywood" ,"cadetblue" ,"chartreuse" ,"chocolate" ,"coral" ,"cornflowerblue" ,"cornsilk" ,"crimson" ,"cyan" ,"darkblue" ,"darkcyan" ,"darkgoldenrod" ,"darkgray" ,"darkgreen" ,"darkkhaki" ,"darkmagenta" ,"darkolivegreen" ,"darkorange" ,"darkorchid" ,"darkred" ,"darksalmon" ,"darkseagreen" ,"darkslateblue" ,"darkslategray" ,"darkturquoise" ,"darkviolet" ,"deeppink" ,"deepskyblue" ,"dimgray" ,"dodgerblue" ,"firebrick" ,"floralwhite" ,"forestgreen" ,"fuchsia" ,"gainsboro" ,"ghostwhite" ,"gold" ,"goldenrod" ,"gray" ,"green" ,"greenyellow" ,"honeydew" ,"hotpink" ,"indianred" ,"indigo" ,"ivory" ,"khaki" ,"lavender" ,"lavenderblush" ,"lawngreen" ,"lemonchiffon" ,"lightblue" ,"lightcoral" ,"lightcyan" ,"lightgoldenrodyellow" ,"lightgreen" ,"lightgrey" ,"lightpink" ,"lightsalmon" ,"lightseagreen" ,"lightskyblue" ,"lightslategray" ,"lightsteelblue" ,"lightyellow" ,"lime" ,"limegreen" ,"line","magenta" ,"mediumaquamarine" ,"mediumblue" ,"mediumorchid" ,"mediumpurple" ,"mediumseagreen" ,"mediumslateblue" ,"mediumspringgreen" ,"mediumturquoise" ,"mediumvioletred" ,"midnightblue" ,"mintcream" ,"mistyrose" ,"moccasin" ,"navajowhite" ,"navy" ,"oldlace" ,"olive" ,"olivedrab" ,"orange" ,"orangered" ,"orchid" ,"palegoldenrod" ,"palegreen" ,"paleturquoise" ,"palevioletred" ,"papayawhip" ,"peachpuff" ,"peru" ,"pink" ,"plum" ,"powderblue" ,"purple" ,"red" ,"rosybrown" ,"royalblue" ,"saddlebrown" ,"salmon" ,"sandybrown" ,"seagreen" ,"seashell" ,"sienna" ,"silver" ,"skyblue" ,"slateblue" ,"slategray" ,"snow" ,"springgreen" ,"steelblue" ,"tan" ,"teal" ,"thistle" ,"tomato" ,"turquoise" ,"violet" ,"wheat" ,"white" ,"whitesmoke" ,"yellow" ,"yellowgreen");

   ########################
   # List of the longer list of colours used by made4 function getcol
   # jmol can't read hexadecimal colours. Tried to convert. But colours not consistent
   ###########################

   @getcols1 =("#E41A1C", "#0000FF","#00FF00", "#FFA500", "#E7298A", "#000080","#00FFFF", "#008000", "#999999", "#800080", "#008080", "#000000", "#FF4500", "#FFFF00","#4DAF4A", "#800000", "#A65628", "#F781BF", "#1B9E77", "#808000", "#377EB8");

   ########################
   # Shorter list of 9 colours used by made4 function getcol
   #######################
   @getcols=("red", "blue", "green", "yellow", "magenta", "cyan", "black","grey", "brown");


   my @cols;
   @cols = @getcols;  # Defines which list from above will be used as colour list

   my @outlist;
   my $i;
   my $no_colours; # number of possible colours
   my $index;
   my $colour;
   my %CHOSEN = ();
 
   # Select $numcol random colours:

   srand;   # initialise random number generator

   $no_colours = $#cols + 1; # number of possible colours
   for ($i = 1; $i <= $numcol; $i++)
   { 
      TRY_AGAIN:
      $index = int rand $no_colours; # gives a number 
                                     # between 0 and ($no_colours - 1)
      if (!($CHOSEN{$index}))
      {
         $colour = $cols[$index];
         @outlist = (@outlist,$colour); # add $colour to @outlist
         $CHOSEN{$index} = 1;
      }
      else
      {   
         goto TRY_AGAIN; # choose another colour as we have this one
                         # already
      }
   }

   return @outlist;
}
  


   


sub display_chime {
  #####################
  # Produce jmol output
  ######################

  my $output_path_base = shift;
  my $classvec = shift;


  my $pdb_filename= $output_path_base. "_li.pdb" ;

  my ($pdbfile) = $pdb_filename =~ /.*\/(.*?)$/;

  `cp $pdb_filename $EP::Config::h->{TMP}/$pdbfile`;

  ## Colour each class with a different color
    my %c;
    $c{$_}++ foreach split /,/, $classvec;
    my $classes = scalar keys %c;

   # my @outcol= javacols($classes);  # javacols not called.

    # Hard-coding in list of colours. Limited to 9 colours 
    # Needs to be fixed.

    my  @outcol=("red", "blue", "green", "yellow", "magenta", "cyan", "black","grey", "brown");
    my @chars = ('A'..'Z');   # Maybe future prob.. limited to 26 groups
    my $i;
    my $index;
    my $str;
    my $str2;
  for ($i=1; $i<=$classes; $i++) {
    $index=$i-1;
    $str = $str."select *".$chars[$index]."; color ".$outcol[$index]."; spacefill 40; ";
    $str2= $str2."Class ".$chars[$index]." is ".$outcol[$index].". ";
  }

  my $embedhtml = <<__EMBED__;
    <script src="$EP::Config::h->{HTML_HREF}/static/js/Jmol.js"></script>
    <script>jmolSetCodebase("$EP::Config::h->{HTML_HREF}/static/js");</script>

     <table cellpadding="0" border="0" width="90%">
     <tr><td>
	  <script>
	    jmolApplet(400, "$EP::Config::h->{HTML_HREF}/tmp/$pdbfile", "wireframe off; select *; spacefill 0; set axes on; color axes black; set ambient 40; background white;set picking ON; set picking LABEL; set perspectivedepth off; set scale3d 1; set frank off; $str");

	  </script>
         <br> "Colours: $str2.  $classvec " <br>
       </td></tr>
       </table>
__EMBED__
  return $embedhtml;
}

1;
